Evaluating the efficiency of coarser to finer resolution multispectral satellites in mapping paddy rice fields using GEE implementation

Timely and accurate estimation of rice-growing areas and forecasting of production can provide crucial information for governments, planners, and decision-makers in formulating policies. While there exists studies focusing on paddy rice mapping, only few have compared multi-scale datasets performance in rice classification. Furthermore, rice mapping of large geographical areas with sufficient accuracy for planning purposes has been a challenge in Pakistan, but recent advancements in Google Earth Engine make it possible to analyze spatial and temporal variations within these areas. The study was carried out over southern Punjab (Pakistan)-a region with 380,400 hectares devoted to rice production in year 2020. Previous studies support the individual capabilities of Sentinel-2, Landsat-8, and Moderate Resolution Imaging Spectroradiometer (MODIS) for paddy rice classification. However, to our knowledge, no study has compared the efficiencies of these three datasets in rice crop classification. Thus, this study primarily focuses on comparing these satellites’ data by estimating their potential in rice crop classification using accuracy assessment methods and area estimation. The overall accuracies were found to be 96% for Sentinel-2, 91.7% for Landsat-8, and 82.6% for MODIS. The F1-Scores for derived rice class were 83.8%, 75.5%, and 65.5% for Sentinel-2, Landsat-8, and MODIS, respectively. The rice estimated area corresponded relatively well with the crop statistics report provided by the Department of Agriculture, Punjab, with a mean percentage difference of less than 20% for Sentinel-2 and MODIS and 33% for Landsat-8. The outcomes of this study highlight three points; (a) Rice mapping accuracy improves with increase in spatial resolution, (b) Sentinel-2 efficiently differentiated individual farm level paddy fields while Landsat-8 was not able to do so, and lastly (c) Increase in rice cultivated area was observed using satellite images compared to the government provided statistics.

backscatter coefficient (measure of reflective strength of radar body) along with Land Surface Temperature (LST) of day and night time and precipitation data can be effectively use to identify the rice fields during the flooding/transplanting and ripening phases with an overall accuracy greater than 80 percent. Data derived from Sentinel-1 provides valuable information in the microwave region and is used to differentiate paddy rice fields in cloud prone regions because of its capability to penetrate cloud cover. On the other hand, areas without cloud cover issue rely on NIR and SWIR data from multispectral instruments. Similarly, GEE platform is delivering users with coarse to moderate resolution imagery, such as Moderate Resolution Imaging Spectroradiometer (MODIS) at 250 m resolution, Landsat-8 at 30 m resolution, and Sentinel-2 at 10 m resolution. These aforementioned satellite instruments are used in previous studies for paddy rice mapping 1,[20][21][22][23][25][26][27] . Sentinel-2 Satellites (2nd series of satellites after Sentinel-1) are part of the Copernicus program initiated by the European Space Agency (ESA) 7 . This satellite mission carries multispectral scanners and consists of two satellites which are Sentinel-2A and Sentinel-2B 28 . In Sentinel-2 satellites, multispectral imaging instruments (MSI) are installed with the ability to record 13 wide-swath bands 29 . Depending upon the band, each satellite provides 10-60 m of spatial resolution with a temporal resolution of 5 days using both satellite constellations. NASA's Landsat-8 mission is a new generation satellite that carries the legacy of Landsat missions. This satellite gets spectral information in visible (V), NIR, SWIR, and Thermal Infrared (TIR) regions 1 . In this satellite mission, two sensors are installed: Operational Land Imager (OLI) and TIR. The OLI obtains spectral information in nine different spectral bands, whereas TIR records thermal information in two bands 30,31 . MODIS mission contains two satellites: Terra and Aqua. These two satellites have nearly the same sensors. The mission monitors the Earth surface every one to two days and obtains spatial data in 36 different spectral bands between 0.405 and 14.385 mm. This mission provides 250 m to 1 km spatial resolution coverage in visible and infrared regions with a 2330 m wide swath 32 .
For preparing paddy rice seasonal spectral profile, Sentinel-2 derived NDVI at 10 m resolution 33 and Sentinel-1 Backscatter coefficient time-series variation 6 were accessed and four observations per month (depending upon satellite revisit intervals) were obtained from 10 randomly selected rice fields inside the Southern Punjab boundary. These training points were then digitized inside GEE, and then their monthly average was calculated individually. For average monthly temperature, MODIS 250 m resolution LST dataset was used 34 , and separate day and night LST were calculated. For precipitation, Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) pentad (5-day) precipitation dataset 35 was used to calculate average monthly precipitation. The average of LST and precipitation were calculated as average values for the southern Punjab inside GEE.  www.nature.com/scientificreports/ Methods. The methodology is divided into three parts. In the first part, rice growth profiles are evaluated using multiple instruments (Sentinel-1, Sentinel-2, Landsat-8, and MODIS), which are crucial to crosscheck training points (non-rice, rice, and water). In second part of methodology, training points (rice, and non-rice) are use to calibrate RF classifier and classify the land cover. After classification, the last part compares derived paddy rice maps using accuracy assessments and differences in the estimated actual areas. The overall methodology performed is shown in Fig. 3. Moreover, all the steps performed in this methodology (i.e., pre-processing and classification) are according to the guidelines and regulations provided by the United States Geological Survey (USGS) (https:// www. usgs. gov), and GEE (https:// www. explo rer. earth engine. google. com/ terms).

Random Forest (RF) based Classification.
Random Forest (RF), a machine learning classifier, has widely been used in several studies for paddy rice mapping [36][37][38][39][40] . Because of its high precision, RF (having a usage percentage of 49% globally in satellite image classification) is the most popular classifier, along with Classification and Regression Tree (13%) and Support Vector Machine (11%) 41 . The RF classifier contains a mixture of tree classifiers in which every individual classifier is created using an arbitrary vector tested independently from the user given input, and each of the trees in the classifier gives a vote for the most excellent prevalent class to classify the given input vector 20 . One of the main advantages of an RF classifier is that it needs two types of parameters to be set, while in contrast, the Support Vector Machine need more than two input parameters by the user. Also, the RF classifier can deal with the data without values, which is not conceivable with Support Vector Machine. Another essential advantage of this classifier is that it also offers the comparative importance of different features during the classification process, which is very valuable in selecting features 42 .
Hyperparameter optimization was performed initially which chooses the optimal parameters to be used for algorithm to yield highest accuracy in a given situation. As a result of hyperparameter optimization, max number of trees value was set to 115 which yield maximum accuracy for our study area. Then the trained classifier was used further to classify paddy rice pixels. For classification validation, the training classifier was fed with validation data, and then the resulting trained classifier was used to calculate the error matrix for resultant maps. All these steps, including cloud masking, data filtering, median, clipping, mosaicking, classifier training, classification, and classification validation, were performed individually for each satellite dataset inside GEE.

Figure 2.
Ground truth data for rice production in southern Punjab for year 2020 (Reported by Crop Reported Service, Government of Punjab). The map is designed in ArcGIS Pro V2.9 software available at ESRI website (https:// www. esri. com/ en-us/ arcgis/ produ cts/ arcgis-pro/ overv iew). The data used in this map is available at CRS website (http:// www. crs. agrip unjab. gov. pk). To map the temporal variation of rice, we used a literature based methodology 6 . The monthly time series of NDVI was created using Sentinel-2 (10 m), Landsat-8 (30 m), and MODIS (250 m) satellite data (Fig. 4). For comparison, the backscatter coefficient (dB) time series was created using Sentinel-1 in dual-band crosspolarization with vertical transmit, horizontal receive (VH) (Fig. 4). The NDVI variation was observed for rice using three multispectral instruments (Fig. 4). Given the rice crop sowing and transplanting phase, NDVI values of Sentinel-2 and Landsat-8 are low (~ 0.16 to 0.20) in months of June and July. However, the values are comparatively higher (~ 0.51 to 0.55) in September and October, representing higher vegetation. For MODIS, peak values of NDVI were observed in August. From Fig. 4, it was observed that backscatter values were low (~ − 18 to − 21 dB) for the flooding stage of rice due to the presence of water and for the harvest stage due to surface scatter from the open land. In contrast, they were significantly high (~ − 16.51 dB) for the vegetation stage due to volume scattering in the rice plant.

Scientific
NDVI variation of various Kharif season crops (Rice, Cotton, and Sugercane) using Sentinel-2 MSI data were compared (Fig. 5a). From their variation, it was observed that NDVI values of crops differ spatially and temporally, therefore, provide great potential in differentiating the crop's growth stages. This phenomenon was  www.nature.com/scientificreports/ also confirmed by previous studies [46][47][48][49][50] . In the last, Fig. 5b shows average day and night LST derived from MODIS satellite along with CHIRPS precipitation data. We used Landsat-8 (100 m thermal data) to derive LST for both day and night. Precipitation time series was prepared using CHIRPS data. CHIRPS is a precipitation dataset which provides 35 years of grided precipitation data with a spatial resolution of 5.5 km 51 . The purpose of LST and precipitation data ( Fig. 5b) was solely to verify the climatic condition for the year 2020 which helped us in planning survey and collecting field samples, in appropriate growth cycle of rice crop in souther Punjab. Figure 5b shows that mean precipitation was high during July to September with moderate mean LST. Literature shows that such weather condition is favorable for enhanced rice growth 52 .
Training samples. Sentinel-1 derived Synthetic Aperture Randar (SAR) data was used to check training points before feeding them to classifier for classification. Previous studies have reported, that the SAR backscatter value is sensitive to water 21 . As paddy rice follows different growth stages, SAR based backscatter can be used to track any pixel either paddy or non paddy. Same principle was applied in this study to differentiate non-rice training points from paddy rice. For specifically non-rice class, each training point was cross checked using Sentinel-1 SAR backscatter data. Only samples which does not follow spectral trend as shown by paddy rice (Fig. 4) were used for representing non-rice training class. To ensure proper coverage and correct selection of the ground truth data for our large study area, we collected the ground truth data (using random sampling approach) through multiple sources: (1) Very High Resolution VHR images from the Google Earth, (2) 167 field rice latitude longitude coordinates point samples randomly collected during the rice field survey in October-November, and (3) Visual interpretation of Sentinel-2 false color composite (FCC) at 10 m spatial resolution. Firstly, we used the Sentinel-2 FCC and VHR Google Earth images to digitize points by visual interpretation except for rice (Fig. 6a). Secondly for each class, we made areas of interest AOIs as circle buffers of the points with the radius of 20 m for Sentinel-2 and Landsat-8 ( Fig. 6b-d). Thirdly, for rice, we took accurate field coordinates  For collecting rice coordinate samples, we choose 6 random locations in the study area which are shown in Fig. 6a (note the 6 areas where rice points are clustered) and collected samples using random sampling approach. For MODIS, however, it was observed that due to its coarser spatial resolution of 250 m, 20 m buffers was not suitable for training points. So, for MODIS, training points without buffer were used and each pixel contributing to paddy rice class was digitized based on its centroid. The AOIs without clear land cover information were omitted. The VHR Google Earth imagery and false-color composite combination for Sentinel-2 at 10 m and field photographs are used to clarify the land cover types, i.e., water, barren, forests, non-paddy, and urban areas. The classes, including non-paddy agriculture, forests, grasslands, urban and barren soil, were merged in one class named non-rice. Water class was kept separate due to its distinct spectral characteristic and was merged with non-rice class after training classifier. The reason for separating water class initially from other land-use classes was to restrict water pixels to create confusion during    53 . In this study, separating water pixels from paddy rice was crucial as they both share similar spectral characteristics 6 (only when rice is in flooding growth stage). The field sizes, shapes, and proximities were considered in labeling land cover information of AOIs. A total of 690 AOIs were collected, including 167 rice samples, 88 water samples, and 435 non-rice samples from different locations of Southern Punjab. The AOIs were distributed randomly, covering the entire study area. The generated AOIs were randomly split (by category) into a 70/30 ratio; 70% were used as training data, while the remaining 30% were used for validation. The paddy rice classification was performed individually for each multispectral satellite dataset utilizing the same training AOIs and the same year (2020).
Pre-processing. Sentinel-2, Landsat-8, and MODIS datasets (freely available inside the GEE data catalog) were used separately and individually for paddy rice classification. Parameters that were kept the same for a fair comparison among these three datasets were (1) same training, and validation data, (2) same satellite data acquisition dates were used (July-October), and the median was taken to ensure best pixel coverage 6 , (iii) Similar spectral bands (if present) were utilized, (4) The same RF classifier was used in each classification, and in the last, (5) same procedure was applied for accuracy assessment and area estimation.
These datasets were used inside GEE. First, these datasets were used with a cloud mask individually, ensuring minimum cloud coverage in data 54 . Then images between July and October were filtered, and their median was taken. Then the datasets were clipped to the boundaries of the study area. For Sentinel-2 and Landsat-8, Visible band with NIR and SWIR bands are used, because NIR and SWIR bands of Sentinel-2 and Landsat-8 enhance crop classification 55 . However, for MODIS, band-1 SWIR and band-2 NIR images with a spatial resolution of 250 m are used 56 .

Results
Rice maps using Sentinel-2, Landsat-8 and MODIS. The rice classified maps for Southern Punjab are generated using Sentinel-2, Landsat-8, and MODIS satellite data (Fig. 8). From their visual comparison, it is observed that MODIS derived maps form clusters instead of classifying individual rice fields, which can be attributed to a coarse-resolution of 250 m. As a result, most pixels of the MODIS map are misclassified because it mixed neighboring crops, and subsequently formed large clusters indicating them rice. Another problem observed in MODIS classified map is denser rice classified pixels clustered in the northeast of southern Punjab, especially in Bahawalnagar, Bahawalpur, and Rahim Yar Khan districts (Fig. 8). This is because of dense vegetation in those areas that can be visually observed in Fig. 1. When compared with Fig. 1, it is confirmed (through visual interpretation) that those areas were actually other crops (Fig. 1) which MODIS has classified totally as rice (Fig. 8).
Landsat-8, on the other hand, performed better than MODIS. At some points, it was successful in classifying individual rice fields. Paddy crop is mostly grown in small, segmented fields; therefore, spatial resolution plays a great role in identifying individual paddy fields 57 . Landsat-8 misclassified some pixels as rice clusters in Khanewal and Dera Gazi Khan districts, that can be clearly seen in the Fig. 8. The rest of the Landsat-8 data is reliable enough for further use. Sentinel-2 classified individual paddy rice with greatest precision as compared to the other two datasets, and successfully differentiating individual paddy rice fields at a resolution of 10 m. The www.nature.com/scientificreports/ performance of these datasets in differentiating paddy rice fields is presented in Fig. 9. For visual comparison, VHR Google Earth imagery (Fig. 9a) and Sentinel-2 derived NDVI at 10 m resolution (Fig. 9b) were used as reference images while derived rice classified maps for Sentinel-2 (Fig. 9c), Landsat-8 (Fig. 9d), and MODIS (Fig. 9e) were compared. At a scale of 1:114,400, a random point (with coordinates 71.37309, 29.94164) was selected and the corresponding derived maps (Fig. 9c-e) were compared at this view. It shows that MODIS formed a large cluster and is not able to identify individual paddy fields. On the other hand, Landsat-8 has identified some individual paddy fields, however, it has misclassified the majority of neighboring small pixels as rice. Comparatively, Sentinel-2 has efficiently classified individual rice fields with lowest number of misclassified pixels. Because of its high resolution of 10 m, Sentinel-2 classified individual paddy fields better than both the Landsat-8 and MODIS.
Technical validation. Through visual inspection (Figs. 8,9) and accuracy assessment (Table 2), we find that the MODIS-derived map is not reliable for differentiating individual paddy fields and cannot be used in  Table 2. Error matrix of paddy rice classification for three satellite data (UA user's accuracy, PA producer's accuracy, OA overall accuracy).

Satellites Class
Error matrix www.nature.com/scientificreports/ local scale regions decision-making. As for the nearly similar area of Sentinel-2 and MODIS datasets, Sentinel-2 derived map was justified through accuracy assessment and visual interpretation. However, MODIS failed to withstand accuracy assessment and visual interpretation mainly due to cluster formation and huge misclassification. Acuracy assessment is important to check the reliability of land use classification and error matrix is a common method to assess the accuracy 58 . An error matrix is also used to draw other statistical measures of accuracy such as overall accuracy (OA), commission error (i.e., user's accuracy (UA)), omission error (i.e., producer's accuracy (PA)), and the kappa coefficient (K) 59 . The UA evaluated by dividing the number of correctly classified pixels (each cateogory) by the total number of classified pixels in that land use category (the row total) 60 . The UA actually shows the probability that the particular pixel of land-use class of classified map actually represent that land-use class on ground. The equation used to derive UA is given as Eq. 1.

UA (%) PA (%) F1-score (%) Kappa (%) OA (%) Non-rice Rice
The PA is computed by dividing the number of correctly classified pixels in each land-use class by the number of reference pixels in that land use category (the column total) 60 . The PA represents how accurately reference pixels of land-use class are classified and the equation used to derive PA is given as Eq. 2.
The overall accuracy (OA) asseses the match or mismatch between classified and ground truth of land use, and is derived by dividing total number of accurately classified pixels by total number of reference pixels 50 . The equation used to derived OA is given as Eq. 3.
The Kappa coefficient evaluates the reliability of classified raster which is calculated using error matrix. The coefficient values vary between 0 and 1, where a higher value means a better reliability and vice versa. The Kappa coefficient is calculated using Eq. 4 50 .
F1-score or F-measure is a useful statistical approach to find the accuracy of classification using precision and recall 61 . F1-score is a harmonic mean of precision and recall, given in Eq. 5.
where precision, also known as the positive predictive value, is calculated by dividing the number of true positive values by the number of all (true and false) positive values (Eq. 6). Recall, also known as sensitivity, is calculated (Eq. 7) by dividing the number of true positive values by total the number of predicted values (true positives plus false negatives) 62 .
The classification process is only considered reliable if it meets accuracy checks, because land use land cover (LULC) maps derived from satellite images may contain some errors due to several factors ranging from techniques in classification to satellite based data retrieval methods 63,64 . For paddy rice maps (Fig. 8), a confusion matrix was developed containing kappa coefficient (K), F1-Score, UA, PA, and OA. It shows that OA was more than 90% for Sentinel-2 and Landsat-8 datasets which supports the reliability of derived maps ( Table 2). For MODIS, OA was 82% showing least reliability for rice mapping compared to Sentinel-2 and Landsat-8.
Area estimation. The area of rice crop is calculated from the classified maps using the "ee.Image.pixelArea()" function available inside GEE. This function computes and assign area for each pixel, and assign area values to each pixel. Then, we apply ee.Reducer.sum which computed total area per class of each classified rice map. For area comparison, crop reported area provided by Crop Reporting Service, Government of Punjab (http:// www. crs. agrip unjab. gov. pk/ repor ts) was used as ground truth data to calculate mean difference between estimated and reported areas (Fig. 10). Multispectral instruments and ground truth data show a good agreement for estimated rice production. The mean percentage difference was approximately 17% for Sentinel-2 and MODIS while 33% for Landsat-8. Previous studies also reported overestimation 37,65,66 , which is mainly because of misclassified pixels in the rice class that may have contributed to some rise in the classified area.
For Landsat-8, the mean percentage difference was more significant (greater than 30%). This overestimation was due to higher misclassified rice pixels as observed previously through visual interpretation of Fig. 8.  www.nature.com/scientificreports/ Surprisingly, the MODIS estimated area was nearly equal to the Sentinel-2 derived area; which reflects why MODIS is used by government officials and many researchers [67][68][69] for large-scale LULC area estimation. To further differentiate the performance of each instrument, and to address the issue why MODIS and Sentinel-2 derived paddy rice area was similar, we compared our study results (fosucing on Multan division) with recently published study performed by Sajjad et al. 70 on crop classification in Multan. The Fig. 11 shows the comparison of paddy rice area from our study for Multan division with area published by Sajjad et al. 70 and Punjab crop statistics report for 2020. The Fig. 11 also shows the mean differences between different area sources compared with the base area estimate provided by Punjab crop statistics report for Multan region. From Fig. 11 it was observed that paddy rice area from Sajjad et al. 70 study using Landsat-8 also showed overestimation (+ 145.7%) Figure 10. Comparison of rice classified cultivated area from GEE and Crop statistics report (local agrarian data). The percentage of overestimation (Δ) by GEE classification relative to ground truth data is provided for each dataset. www.nature.com/scientificreports/ whereas our derived area using Landsat-8 also shows similar overestimation (+ 143.5%). Furthermore the Fig. 11 also highlighted a major difference in area estimation of Sentinel-2 (+ 26.3%) and MODIS (+ 165.5%) compared with Punjab crop statistics report.

Discussion
Rice mapping is essential for bringing water-use efficiency, achieving regional food security, controlling disease transmission, and countering global warming. As our study area is a semi-arid region with lower precipitation (Fig. 5b), free available and GEE integerated MSI are suitable to map paddy rice crop due to lower cloud cover in the cropping months 27,71,72 . Despite their frequent use in mapping paddy fields, no study has compared the efficiency of Sentinel-2, Landsat-8 and MODIS datasets in classifying individual paddy rice fields. Our study compares these datasets by estimating their paddy crop classification potential, using accuracy assessment methods, and the mean difference between estimated area and area based on ground truth data.
The NDVI values derived from Sentinel-2 at 10 m spatial resolution (Fig. 5a) were used to track changes in crop growth stages. The NDVI of any crop differs spatially and temporally and thus can be used as a valuable tool to study crop growth changes over space and time 48,49,71,73,74 . The temporal variation shows that rice sowing and transplanting stages start from May until June in southern Punjab region. Lower NDVI values (Fig. 5a) in July and August reflect significant NIR adsorption on the water surfaces which in this case are paddy rice fields flooded with water and with minimum green area. Similarly, Usman et al. 50  The class 'rice' exhibits a unique trend compared to other crops (Cotton, and Sugarcane) in the Kharif season (Fig. 5a). For rice crop, the initial NDVI variation is a bit slow and lengthy because of rice nursery growth in June. However, the latter part of NDVI variation reaches its maximum due to rapid rice crop growth in September and October. Figure 7b shows variation in temperature and precipitation patterns during the Kharif season. It was observed that from July to September, mean precipitation values were high along with moderate LST (Day and Night), which supported rice cultivation 14,52 .
For rice crop, Sentinel-2 and Landsat-8 data showed similar NDVI trends (Fig. 4). Due to their high spatial and temporal resolution, these trends were successful in differentiating various rice growth stages. Conversely, the MODIS curve showed a peak response in August, mainly due to its coarse resolution (of 250 m). Although MODIS has a higher revisit time (1 day) than Sentinel-2 (5 days) and Landsat-8 (16 days), it cannot be used to study the growth cycle of individual rice fields due to its coarser resolution. Likewise Chen et al. 75 concluded that the MODIS imagery with 250 m resolution could not provide accurate vegetation information on corn growth. The variation in Sentinel-1 backscatter coefficient for paddy rice (Fig. 4) showed that the coefficient is very sensitive to water in paddy rice fields. The backscatter values were low in June and July due to the presence of water in the fields, and were high in September and October due to high vegetation. Singha et al. 6 noted a similar trend in the the backscatter values. The Sentinel-1 VH backscatter coefficient can be used as an indicator for tracking paddy rice growth as the backscatter values change with the varying conditions of paddy rice stages. Unlike other crops, rice spends a substantial time submerged in water. Thus, it can be easily differentiated from other crops. Also, it has distinct growth stages, which can be observed using the Sentinel-1 VH backscatter coefficient. Based on the temporal variation of NDVI and Sentinel-1 backscatter coefficient, an accurate rice growth profile was generated for the Southern Punjab region for the year 2020, which was used as a baseline in rice fields sampling. The coordinates used for rice training and validation by global positioning system (GPS) were crosschecked using NDVI and backscatter temporal variation. This step was performed inside GEE to ensure that only rice fields coordinates having the same trend as observed in Fig. 5 were taken for further classification and accuracy assessment. Previous studies 6,20,76-78 support the RF classifier's capability in identifying paddy fields. Therefore, this study utilized RF classifier for the classification of paddy rice.
Results showed that Sentinel-2 performed best in classifying individual rice fields (Figs. 8,9). MODIS failed to accurately classify rice fields, and it formed clusters with neighboring fields. A visual interpretation was performed for comparing the efficiency of each dataset in classifying individual rice fields. It was observed that Landsat-8 efficiently differentiated paddy rice fields; however, it misclassified rice pixels at some points with other crops having similar vegetation characteristics. Overall, Sentinel-2 performed best in accuracy assessment than both Landsat-8 and MODIS (Table 2).
Further, we compare estimated and reported areas for rice production, which shows an overestimation of ~ 18% for both Sentinel-2 and MODIS, and ~ 34% overestimation for Landsat-8 (Fig. 10). Similar overestimation was reported by other researchers 37,65,66 . For example, Mananze et al. 37 used GEE to map shifting agriculture dynamics in Mozambique and found 15.5% of mean difference between the estimated and ground truth areas for agricultural land use. In another study to map cropping intensity of small-scale farms in India, Jain et al. 79 reported a large discrepancy (with percentage difference of 150.8%) between cropland area estimated through remote sensing (Landsat data) and the area obtained from agricultural census. Likewise, Dheeravath et al. 80 found that irrigation area derived using MODIS (500 m) data was greater than the ground truth data (with the mean percentage difference of 17.2%) in most Indian States.
To compare the overestimation in paddy rice area, we compared our results with a previous study performed by Sajjad et al. 70 focusing on Multan division (Fig. 11). The Fig. 11 shows area estimated as well as mean percentage difference of area from each source compared with crop statistics report. It was observed that there exist a similarity in overestimation of Landsat-8 derived area reported by Sajjad  The multispectral imagery-based paddy rice mapping implemented in GEE includes several benefits. First, there is reduction in data acquiring and pre-processing time 40 . Otherwise, it might have taken longer to download and pre-process the datasets for a large study area like ours (i.e., southern Punjab). Second, there is a significant decrease in computational time 81 . The MSI-based RF algorithm inside the GEE performed paddy rice fields classification in a very short time. Thus, high-performance computing resources like GEE facilitate quick and rapid mapping of paddy rice planting areas at a large scale 1,6 . Future studies of paddy rice mapping under similar climatic conditions can utilize the most efficient dataset found in this study: Sentinel-2. In addition, further research may focus on rice mapping for the entire country using high-resolution Sentinel-2 imagery with a similar GEE-based methodology.

Conclusions
In this study, we compared three, free and open, satellite datasets (Sentinel-2 (10 m), Landsat-8 (30 m), and MODIS (250 m)) to evaluate their efficiency/effectivity for mapping of rice crop in southern Punjab, Pakistan. The results show that rice classification improves with increase in spatial resolution. Among the three, the coarsest resolution dataset, MODIS, do not precisely classify individual rice fields. However, the other two highresolution datasets (i.e., Sentinel-2 and Landsat-8) can accurately track crop growth stages for making its growth profiles and studying temporal trends. From the visual interpretation, accuracy assessment, and mean difference in area estimation, we observe that MODIS was not able to classify individual rice fields correctly in this case. However, Landsat-8 classified the paddy rice fields accurately with moderate misclassification. In comparison to MODIS and Landsat-8, Sentinel-2 performed best in differentiating individual rice fields with lower misclassification and showed a lesser mean difference in area estimation. Due to limited resources and training data, this study focused only on the RF classifier to classify southern Punjab divisions. JavaScript API is used inside GEE, which facilitates future studies to implement the same methodology in GEE with ease. Future studies may focus on (1) rice mapping of the whole of Pakistan, which could be easily implemented inside GEE using the same methodology and code, (2) comparison of different machine learning classifiers (like Random Forest RF, Support Vector Machine, and Artificial Neural Network to classify paddy rice efficiently, and (3) comparison of Synthetic Aperture Radar and Multispectral Instruments to map paddy rice efficiently.

Data availability
All the datasets used in this study are open access, and references are provided at first mention. Paddy rice classification images for Southern Punjab are also available, and can be provided upon reasonable request to corresponding authors. No direct plant material is used in this study. Only open access multispectral satellite images are used, that are freely available inside GEE data catalogue (https:// www. devel opers. google. com/ earthengine/ datas ets).